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ABSTRACT 

This  paper  considers  the  problem  of  tracking  and  intercepting  a  potentially  moving  ground-based  RF  source  with  an 
autonomous  guided  munition  that  has  a  passive  bearings-only  sensor  located  on  its  nose.  It  is  assumed  that  the  munition 
has  lost  GPS  signal  lock  and  that  it  relies  only  on  its  noisy  inertial  measurement  unit  (IMU)  for  guidance  and 
navigation.  Bearings-only  target  motion  analysis  (TMA)  algorithms  are  used  to  obtain  a  position  and  velocity  estimate 
for  the  RF  source  using  the  position,  velocity  and  attitude  estimates  of  the  munition  as  well  as  the  azimuth  and  elevation 
measurements  obtained  from  the  bearings-only  RF  sensor.  Six  degree-of-freedom  (6-DOF)  and  three  degree-of-freedom 
(3-DOF)  munition  models  are  used  to  evaluate  the  tracking  and  intercept/seeker  performance  of  a  hybrid  coordinate 
(HC)  extended  Kalman  filter  (EKF),  a  particle  filter  (PF),  a  multiple  hypothesis  (MF1)  FIC-EKF,  and  a  pseudo-linear 
least  squares  (PLLS)  filter. 
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1.  INTRODUCTION 

Civilian  and  military  navigation  systems  are  increasingly  relying  on  the  accuracy  of  the  information  provided  by  the 
Global  Positioning  System  (GPS).  Moreover,  since  GPS  receivers  are  now  embedded  in  many  different  types  of 
systems,  this  dependence  goes  far  beyond  inertial  navigation  and  guidance  systems  and  extends  into  areas  such  as 
personal  communication  systems  (PCS)  and  wireless  internet  access  systems.  In  military  systems,  the  loss  of  GPS  signal 
lock  could  cause  an  otherwise  successful  mission  to  become  a  failure,  thus  wasting  valuable  resources.  In  some  cases,  it 
could  even  mean  the  difference  between  life  and  death.  The  same  is  true  in  civilian  applications.  Therefore,  it  becomes 
imperative  not  only  to  protect  the  integrity  of  the  GPS  signal,  but  to  locate  and  to  eliminate  any  threats  to  GPS  as  soon 
as  possible. 

Although  GPS  has  some  built-in  anti-jam  (AJ)  capabilities,  it  can  be  further  augmented  by  integrating  a  GPS  receiver 
with  an  Inertial  Navigation  System  (INS)  and  by  utilizing  AJ  GPS  antennas.  The  performance  of  GPS/INS  systems  is 
especially  crucial  for  guided  missiles  with  extended  ranges.  Depending  on  the  variant,  the  standoff  range  of  unpowered 
guided  missile  systems  can  be  as  much  as  100  kilometers.  When  a  guided  missile  is  tasked  to  fly  deep  behind  enemy 
lines  one  should  expect  jammers  to  be  placed  anywhere  along  its  trajectory,  thereby  increasing  the  probability  of  losing 
GPS  signal  lock  sometime  during  the  mission.  If  the  flight  time  of  an  unpowered  weapon  is  more  than  1 00  seconds  it  is 
in  danger  of  not  achieving  its  accuracy  requirements  due  to  the  potential  loss  of  GPS  signal  lock  somewhere  along  its 
trajectory.  This  is  because  once  GPS  signal  lock  has  been  lost,  navigation  errors  grow  at  a  rate  which  depends  on  the 
quality  of  the  INS.  A  one  degree-per-hour  INS  will  have  about  a  30  meter  circular  error  probable  (CEP),  that  is,  the 
radius  which  encloses  50%  of  all  samples,  within  100  seconds  [1,  2].  This  situation  creates  a  dilemma  for  mission 
planners.  That  is,  what  is  the  weapon  to  do  if  it  has  lost  GPS  lock  and  is  still  too  far  away  from  its  target  to  meet  its 
accuracy  requirements?  One  solution  is  to  equip  the  munition  with  an  angle-of-arrival  (AO A)  sensor  that  can  provide 
the  mission  planner  with  a  number  of  new  opportunities  in  the  event  that  GPS  signal  lock  is  lost.  First,  if  so  desired,  it 
provides  the  missile  with  a  jammer  homing  system.  This  is  especially  useful  if  the  missile  determines  that  the  jammer 
and  the  designated  target  are  collocated.  In  which  case,  it  will  permit  the  weapon  to  successfully  complete  its  mission 
with  potentially  improved  accuracy.  Second,  by  sharing  this  information  with  other  cooperative  missiles  and/or  a  central 
command,  it  will  improve  the  success  probability  of  future  follow-on  missions.  Finally,  the  angle-of-arrival  and  power 


*  kezal@toyon.com;  phone  1  805  968  6787  xl80;  fax  1  805  685  8089;  www.toyon.com 


Report  Documentation  Page 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1.  REPORT  DATE 

2006 

2.  REPORT  TYPE 

3.  DATES  COVERED 

00-00-2006  to  00-00-2006 

4.  TITLE  AND  SUBTITLE 

5a.  CONTRACT  NUMBER 

i racking  ana  interception  ot  grouna-nasea  Kr  sources  using  autonomous 
guided  munitions  with  passive  bearings-only  sensors  and  tracking 

5b.  GRANT  NUMBER 

algorithms 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Toyon  Research  Corporation, 75  Aero  Camino  Suite  A,Goleta,CA,93117 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

The  original  document  contains  color  images. 

14.  ABSTRACT 

see  report 

15.  SUBJECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 
ABSTRACT 

18.  NUMBER 
OF  PAGES 

12 

19a.  NAME  OF 
RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


level  information  obtained  from  multiple  missiles  within  a  common  area-of-interest 
can  be  fused  together  to  improve  our  battlefield  and  situational  awareness.  With  this 
added  knowledge  it  will  be  possible  to  compute  trajectories  which  could  take 
advantage  of  regions  of  lower  jamming  power  levels  to  reach  designated  targets. 
Power  levels  recorded  by  prior  missions  can  also  be  used  to  determine  the  locations 
of  large  jamming  sources,  thus  providing  the  location  of  new  potential  targets. 

We  are  primarily  interested  in  evaluating  the  performance  of  bearings-only  target 
motion  analysis  (TMA)  algorithms  that  use  munition-based  AOA  sensor 
measurements  as  a  means  of  tracking  and  homing  on  GPS  jammers.  The  bearings- 
only  tracking  problem  is  also  known  as  a  passive  tracking  problem  due  to  the  fact 
that  active  radar  is  not  utilized  to  obtain  measurements  of  the  target’s  motion. 
Typically,  the  measurements  are  the  bearing  angles  in  the  direction  of  the  target 
emissions  in  the  form  of  elevation  (9)  and  azimuth  (4>)  angles.  The  emissions  may  be 
any  form  of  RF  emissions  such  as  communication  signals  and  radar  or,  perhaps, 
sonar  signals  emanating  from  the  target. 


Figure  1:  Scenario 
geometry 


Bearings-only  TMA  is  a  difficult  problem  to  address  because  it  has  inherent 
observability  issues  [3,  4,  5]  and  the  system  equations  have  nonlinear  components 
making  it  difficult  to  construct  stable  and  unbiased  estimators.  In  the 
standard  extended  Kalman  filter  (EKF)  approach  to  this  nonlinear 
problem  the  target  state  is  defined  in  Cartesian  coordinates  so  that 
the  propagation  of  the  state  (i.e.  the  dynamic  equation)  is  linear, 
while  the  nonlinear  measurement  equation  is  linearized  through  a 
Taylor  series  expansion.  The  Cartesian  coordinate  EKF  (CC-EKF) 
is  known  to  have  divergence  issues  for  the  bearings-only  problem, 
particularly  when  the  target  initialization  is  poor  [3].  A  particular 
variation  of  the  CC-EKF  is  the  Modified-Gain  EKF  (MG-EKF)  of 
Song  and  Speyer  [6]  which,  in  our  implementation,  had  similar 
divergence  issues.  Another  algorithm,  dubbed  the  pseudo-linear 
least-squares  (PLLS),  is  obtained  through  some  manipulation  of  the 
measurement  equation  and  has  very  good  convergence  properties 
for  stationary  targets  (Aidala,  Nardone)  [7].  However,  the  PLLS 
does  lead  to  a  biased  estimator  [8].  In  this  algorithm  the  set  of 
pseudo-linear  measurements  are  linear  in  the  Cartesian  state.  The 
PLLS  measurement  matrix  is  a  function  of  the  range  to  the  target 
and  includes  noise  which  ultimately  leads  to  a  bias  in  the  state  estimates. 


Figure  2:  Coordinate  system 


Another  variation  of  the  EKF  is  the  Modified  Polar  Coordinate  EKF  (MP-EKF)  where  the  state  dynamics  are  nonlinear 
and  the  measurement  equation  is  linear  [9].  The  primary  drawback  to  the  MP-EKF  approach  is  the  complexity  of  the 
nonlinear  state  equation  which  requires  an  approximation  for  the  state  estimate.  In  contrast,  the  hybrid  coordinate  EKF 
(HC-EKF)  is  based  on  the  observation  that  the  state  dynamics  are  linear  in  Cartesian  coordinates  while  the  measurement 
equation  is  linear  in  modified-polar  coordinates  [10,  1 1].  The  approach  is 
then  to  implement  the  Kalman  filter  equations  in  the  coordinate  system 
for  which  the  process  is  linear.  Thus,  for  propagation  of  the  state 
estimate  and  error  covariance  matrix,  the  Kalman  filter  equations  are 
implemented  in  Cartesian  coordinates,  and  for  state  and  covariance 
updating,  the  Kalman  filter  equations  are  implemented  in  modified-polar 
coordinates.  The  transformation  of  the  state  estimate  from  one 
coordinate  system  to  the  other  is  an  exact  mathematical  process;  whereas 
transforming  the  covariance  requires  an  approximation.  We  found  that 
the  HC-EKF  preserves  the  stability  of  the  MP-EKF  formulation  while 
retaining  the  simplicity  of  the  CC-EKF.  Recent  research  has  attempted  to 
improve  the  performance  of  the  CC-EKF  by  considering  multiple  initial 
conditions  using  parallel  Kalman  filters  and  weighting  the  outcome  of 
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Figure  3:  INS  error  model 
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each  filter  to  come  up  with  a  single  estimate  [12,  13].  We  implemented 
the  multiple-hypothesis  approach  with  the  HC-EKF  but  did  not  see  any 
improvement  in  performance.  This  was  mostly  likely  due  to  the  fact 
that  the  HC-EKF  is  inherently  stable  and  has  good  convergence 
properties.  Another  family  of  filtering  algorithms  which  has  gained  - 
much  interest  in  the  last  decade  is  the  group  of  algorithms  referred  to  as  | 
particle  filters  [14].  Particle  filters  are  based  on  Monte  Carlo  simulation 
methods.  One  advantage  of  particle  filters  is  that  there  is  no  assumption  f 
of  linearity  or  Gaussian  noise  processes.  The  probability  density  J 
function  (PDF)  of  the  jammer  is  approximated  with  points  (the 
particles)  in  the  state  space.  As  the  number  of  points  representing  the 
PDF  increases,  the  computational  expense  associated  with  the 
algorithm  increases  as  well. 

In  this  paper  we  consider  a  simple  guided  munition  scenario  with  an 
AOA  sensor  on  the  nose  of  the  munition.  The  scenario  is  described  in 
the  following  section.  We  also  limit  ourselves  to  the  discussion  of  the  relative  performance  of  the  PLLS  algorithm,  the 
HC-EKF,  the  multiple-hypothesis  HC-EKF  (MH-HC-EKF),  and  the  particle  filter  (PF)  implementation.  These 
algorithms  are  briefly  described  in  Section  3.  We  provide  some  typical  simulation  results  in  Section  4  and  concluding 
remarks  in  Section  5. 

2.  SCENARIO  DESCRIPTION  AND  Sensor 

ASSUMPTIONS 

A  bird’s  eye  view  of  the  guided  munition  scenario  is 
shown  in  Figure  1.  A  flat  earth  model  is  used  with  an 
east-north-up  (ENU)  inertial  coordinate  system.  The 
sensor/munition  coordinate  system  is  defined  in 
Figure  2  with  its  x-axis  going  forward  through  the 
nose  of  the  munition.  The  inertial  position  of  the 

guided  munition  is  denoted  by  \M  (t)  e  R  .  The 

intended  target  is  at  the  origin  and  the  jammer  inertial 

..  ..  ,,  i  T  .  .  ,  Figure  5:  Simulation  block  diagram 

position  is  denoted  by  x  j  (t)  e  R  .It  is  assumed  that 

the  jammer  is  traveling  with  a  constant  velocity  \(t)  =  \'  j  .  Our  simulation  uses  two  different  models  for  the  missile 
dynamics.  Both  models  attempt  to  capture  the  missile  dynamics  of  an  extended  range  munition  and  the  navigation  errors 
of  both  models  are  the  same.  The  INS  error  model  is  shown  in  Figure  3.  The  first  missile  model  is  a  3-DOF  dynamic 
model  that  computes  the  axial  and  normal  aerodynamic  forces  acting  on  the  body  of  the  munition  through  a  set  of  look¬ 
up  tables  that  are  a  function  of  the  Mach  number  and  the  angle  of  attack  of  the  weapon.  A  piecewise  exponential 
atmospheric  density  model  is  employed.  Attitude  error  is  modeled  by  adding  rotational  errors  that  are  consistent  with 
the  INS  error  model  of  Figure  3.  The  second  model  is  a  full-blown  6-DOF  dynamic  model  obtained  from  the  Air  Force 
using  the  VisSim™  environment  by  Visual  Solutions  Corporation. 

The  munition  is  released  at  time  t  =  TL  <  0  and  loses  GPS  signal  lock  at  time  t  =  0.  We  are  primarily  interested  in  the 
period  of  time  after  GPS  signal  lock  has  been  lost.  Once  GPS  signal  lock  has  been  lost  the  weapon  can  do  one  of  two 
things:  it  can  continue  to  fly  towards  its  intended  target  or  it  can  home  on  the  jammer.  At  t  —  0  the  munition  begins  to 
utilize  the  AOA  sensor  located  on  the  nose  of  the  weapon.  The  sensor  measures  the  elevation  angle  (0)  and  azimuth 
angle  (()))  of  the  jammer  relative  to  the  munition  in  the  sensor/munition  coordinate  system.  Bearing  measurements  are 
taken  every  10ms  and  these  samples  are  averaged  every  100ms.  The  elevation  angle  dependent  Gaussian  sensor  noise 

model  used  in  our  simulations  is  represented  in  Figure  4  by  the  composite  angle  error  cr“  =  <Jq  +  sin  (6)  cr^  where 
E{uq}  =  =  0  and  E{hq}  =  Ety^}  =  ctq  =(7^  .  The  equation  relating  these  angle  measurements  is  given  by 
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Figure  4:  Angle  error  model 
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where  r(£)ef?  is  the  unit  bearing  vector  in  the  direction  of  the  jammer  defined  in  the  sensor/munition  coordinate 
system. 

The  bearing  measurements  are  used  to  compute  an  estimate  of  the  jammer  location  using  bearings-only  localization 
algorithms.  Selected  algorithms  are  described  in  Section  3.  The  estimated  jammer  location  and  our  confidence  in  the 
estimate  are  used  to  determine  whether  or  not  the  munition  will  home  on  the  jammer  or  continue  toward  its  intended 
target.  The  decision  algorithm  itself  is  outside  the  scope  of  this  paper  and  is  not  discussed.  Figure  5  shows  the  block 
diagram  of  the  simulation  components. 
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ALGORITHM  DESCRIPTIONS 
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weapon  in  Cartesian  coordinates.  Hence,  x^  =  —  x = 
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where  a  j  represents  the  acceleration  of  the 


jammer,  and  a0  represents  the  acceleration  of  the  weapon  (i.e.  sensor  platform)  which  is  assumed  to  be  measured  by 
the  INS.  We  can  put  the  above  equation  in  matrix-vector  form  as 
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where  I3  is  the  three  by  three  identity  matrix  and  is  the  three  by  three  zero  matrix.  Since  the  propagation  is  done 
in  discrete  time  steps,  the  above  equation  can  be  discretized  by  integrating  over  one  time  step  (T).  Furthermore,  we  have 
assumed  that  the  jammer  is  not  accelerating,  so  that  a7  ~  0  .  The  resulting  discrete-time  equation  modeling  the  system 
dynamics  is  given  by 
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a 0(k)  or,  equivalently,  xc(k  + 1)  =  ®^xc(A:)  +  r^.ulk) . 


Our  goal  is  to  estimate  the  position  of  the  jammer  through  an  estimate  of  the  range:  Xj  =  x  M  +  R  .  An  estimate  of  the 
missile  position,  xM  ,  is  available  from  the  INS. 


The  measurement  equation  (CCME)  can  be  rewritten  as 


z  c(k) 


(^) 

0m(k) 


=  f(A(Jfc)r(Jfc))  +  n(Jt) 


(CCME) 


3  vT 

where  A(k)  e  R  is  a  rotation  matrix  which  transforms  the  unit  bearing  vector  in  the  inertial  coordinate  system 
r(k )  e  R  to  the  unit  bearing  vector  in  the  antenna  coordinate  system  r  (k) .  Since  f(«)  is  a  nonlinear  function,  and 
r(k)  =  R(£)/|R(£)| ,  the  measurement  equation  (CCME)  is  a  nonlinear  function  of  the  state  x^ .  Hence,  it  is  difficult  to 

obtain  an  estimate  of  the  state  without  some  form  of  approximation.  In  fact,  it  has  been  demonstrated  that  the  extended 
Kalman  filter  in  Cartesian  coordinates,  which  uses  a  Taylor  approximation  of  (CCME),  can  lead  to  instabilities  in  the 
estimation  process  [3],  Selected  alternative  bearings-only  estimation  algorithms  are  described  in  the  following  sections. 


3.1 


Pseudo-Linear  Least  Squares  (PLLS)  Algorithm 


In  its  simplest  incarnation,  the  PLLS  algorithm  assumes  that  the  jammer  is  stationary,  that  is  Xj  =  0  .  Since  the  jammer 
is  assumed  to  be  stationary,  there  are  no  dynamics  to  update  and  the  observations  are  assumed  to  be  pseudo-linear. 
While  this  fdter  is  very  simple  to  implement  and  remains  stable,  it  is  a  biased  estimator  [8],  The  filter  is  derived  by 
defining  bearing  measurements  in  inertial  coordinates.  Px/z(k)  =  arctan(R  x{k)/Rz{k))  =  arctan(rA'(^)/rz(^)) 
with  Py/x(k)  and  Pz/y(k)  defined  in  a  similar  manner  and  where 

XJX  (k)  ~  XMX  W  rx(k) 

x  (k)  ~  XM;,  (k)  =  R Y(k)  =|R(*)|  rY(k)  =  |R(^)|Ar  (/t)r  (/t) . 

*jz(k)-xMz{k)  |_Rz(^)J  \jz(k)_ 

Assuming  an  additive  Gaussian  noise, 

Px/z(k)  arctan(Rx(£)/Rz(£))  n  x/z(k)  arct<m{rx(k)/rz(k)) 

J3v/x(k)  =  arctan(R7(&)/R;^ (£))  +  n  v/x(k)  =  arctan(ry  (£)/rx (A:))  +  n(£) ,  with  cos(n(£))  «  1 , 
Pz/y(k)  |_  arctan(Rz  (&)/Ry  (&))J  nz/y(k)  larc\an{rz(k)lrY(k))_ 

and  defining  Rxlz (k)  =  <Jr]i(k)  +  r %(k)  ,Ry/x(k)  =  Jrx(k)  +  ry(k) ,  Rzly (k)  =  -Jryik)  +  rz(k) , 
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we  can  write  the  measurement  equation  as  a  linear  function  of  the  inertial  state  of  the  jammer  and  the  munition: 

Zp(k)  =  HT(k)xM  =nT(k)xj  -nP(k).  (PLME) 

T 

We  note  that  zP(k)  =  H  (k)x ^  is  a  measured  quantity  since  it  is  a  function  of  the  (noisy)  bearing  angles  in  inertial 
coordinates.  The  term  njo(£)  is  the  effective  noise  of  the  system.  The  measurement  equation  (PLME)  can  now  be 
directly  applied  to  the  standard  least-squares  recursive  form: 

K(k)  =  P (k  -  l)H(/t)[l  +  Hr  (k)P(k  -  1)H(£) 

x j(k)  =  x j {k  -\)  +  R{k'}^. P{k) -RT {k)x j(k- 1)  (PLLS) 

P(Jfc)  =  P(Jfc  - 1)  -  K(yt)Hr  (Jfc)P(*  - 1) 

where  x  j  (0)  =  Xq  is  the  initial  estimate  of  the  jammer  location  and  P(0)  =  21  >  0  for  any  positive  constant  'k. 

In  our  simulations,  the  covariance  P  is  reset  every  5  seconds  with  X  =  0.02,  and  x0  is  initialized  after  4  seconds  of 
filtered  bearing  measurements  have  been  accumulated.  The  filtered  bearing  measurements  are  extrapolated  to  intersect 
ground  level  (defined  as  the  target  altitude)  to  determine  the  initial  position  estimate.  The  separation  angle  between  the 
latest  LOS  vector  to  the  jammer  position  estimate  is  continuously  monitored.  The  filter  is  reinitialized  if  the  separation 
angle  exceeds  five  degrees.  The  new  initial  condition  is  selected  along  the  latest  measured  LOS  vector  at  a  range  equal 
to  the  last  position  estimate  range. 

3.2  Hybrid  Coordinate  Extended  Kalman  Filter  (HC-EKF) 

The  HC-EKF  is  based  on  the  observation  that  the  measurement  equation  is  linear  in  modified-polar  coordinates  while 
the  state  dynamics  are  linear  in  Cartesian  coordinates.  The  approach  is  then  to  implement  the  Kalman  filter  equations  in 


the  coordinate  system  for  which  the  process  is  linear.  The  jammer  state  in  modified-polar  coordinates 
1  “iT 

\Mp  =  |R|  ( (c/|R|  / i/?)/|R|)  e  R 8  includes  the  unit  line-of-sight  (LOS)  vector  xMi>  =  r  sR 3  ,  the  time- 

derivative  of  the  unit  LOS  vector  x MPi  =  r  eR  ,  the  inverse  of  the  range  x MP}  =  |R|  e R  ,  and  the  range  rate  divided 
by  the  range  \mp4  =  (c/|R|  /  cfr)/|R|  eR .  In  the  standard  hybrid-coordinate  EKF  implementation,  the  measurement  is  a 

direction  cosine  vector  in  the  direction  of  the  jammer  and,  hence,  the  measurement  equation  in  modified-polar 
coordinates  is  a  linear  function  of  the  state  [10,  11]: 

zMp(k)  =  Aix3(k)r(k)  +  n(k)  =  [\ixi(k)  03x3  0  o]xMP(k)  +  n(k) .  (MPME) 

Thus,  for  propagation  of  the  state  estimate  and  error  covariance  matrix,  the  Kalman  filter  equations  are  implemented  in 
Cartesian  coordinates,  and  for  state  and  covariance  updating,  the  Kalman  filter  equations  are  implemented  in  modified- 
polar  coordinates.  The  transformation  of  the  state  estimate  from  one  coordinate  system  to  the  other  is  an  exact 
mathematical  process  whereas  transforming  the  covariance  requires  a  Jacobian  and,  hence,  an  approximation. 

In  our  implementation,  we  note  that  a  unit  vector  in  the  direction  of  the  measured  LOS  is  equivalent  to  measuring  the 
elevation  and  azimuth  angles.  If  the  measurement  were  the  bearing  vector  and  not  the  bearing  angles,  the  resulting 
measurement  equation  would  be  linear  in  modified-polar  coordinates.  In  fact,  the  measurement  would  be  a  noisy 
version  of  xMPi  ( k ) .  Thus,  we  transform  the  measurements  from  bearing  angles  to  a  bearing  vector  representing  a  unit 

LOS  vector  to  the  jammer: 

cos(dm(k)) 

zH(k)  =  sin(dm(k))cos((/>m(k))  *  [A(£)  :  03x5]xMp  +  n(£)  =  r(k)  +  n(£)  =  A(k)r(k)  +  n(£)  (HCME) 
sm(0m(k))sin(<fm(k))  _ 

where  03x5  is  the  3  by  5  zero  matrix.  Note  that  the  equation  above  is  an  approximation  since  the  noisy  measured 
elevation,  9m  ( k ) ,  is  given  by  6m  ( k )  =  0{k)  +  np(k)  where  6(k)  is  the  true  elevation  angle.  In  order  to  use  the  Kalman 
filter  update  equation  it  is  necessary  to  derive  the  approximate  noise  characteristics  of  n (k) ,  or  the  measurement  error 
covariance  matrix.  Hence,  we  must  expand  the  sinusoidal  functions  and  find  the  mean  and  variance  of  the  “noise”  term 
for  each  component  of  the  measurement  vector.  For  example,  we  expand  the  cosine  function  for  the  first  component  of 
zH(k)  as  follows 

z\(k)  =  cos(9(k))cos(np(k))  -sin(9(k))sin(/7P(k)) 

2 

where  ng(k)  is  a  zero-mean  Gaussian  random  variable  with  a  variance  of  Oq  .  Assuming  that  the  magnitude  of  the 
noise  is  small,  we  approximate  the  cosine  of  the  noise  as  one  which  yields 

zy(k) «  cos(<9(£))+  n\(k) 

where  n\{k)  =  sin(<9(£))sin(«0(A:)) .  It  can  be  shown  that  the  mean  and  variance  of  n\{k)  are  given  by 

E{nx{k)}  =  Q 

E[il(k)}=  a[(k)  =  3-sin2  (d{k))  1  -  e~2f7 » 

In  the  same  way,  we  can  determine  the  noise  statistics  for  z^{k)  and  z3(£)  : 

z2(k)  =  sin(#OT (£))cos(^m (£)) «  sin(0(£))cos(^(£))+  n2(k) 
z3 (k)  =  sin (0m (£))sin(^m (k)) «  sin(^(A:))sin(^(A:))+  n3(k) 


The  statistics  on  n2(k)  and  n^k)  are  assumed  Gaussian  with 
E{n2(k)\=  E{n3(k)\=0 


Em 2  {k)\=  a\ (k)  =  0.25cos2 (0{kj) cos2 {<f(k))  1  - e~2aS 
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We  stress  that  the  actual  measurements  are  the  elevation  and  azimuth  angles  representing  the  direction  from  which  the 
jammer  energy  emanates.  These  measurements  are  transformed  so  that  the  measurement  equation  can  be  expressed  as  a 
linear  function  of  the  modified-polar  state  variables.  Note  that  the  variance  on  the  angle  estimates  is  a  function  of  the 
actual  angles.  For  example,  the  variance  on  the  elevation  angle  estimate  increases  as  the  elevation  angle  increases.  Thus, 
throughout  the  flight  of  the  weapon,  the  accuracy  of  the  angle  of  arrival  estimates  varies,  and  we  want  our  measurement 
noise  covariance  to  reflect  this  difference.  Thus,  the  measurements  are  weighted  according  to  their  accuracy.  By 
accounting  for  the  true  time -varying  nature  of  the  measurement  error  covariance  matrix,  we  build  a  more  robust 
algorithm. 

The  transformations  between  the  Cartesian  and  modified-polar  coordinate  systems  are  derived  in  a  straightforward 
manner;  therefore,  only  the  results  will  be  given  here  (see  [10  and  11]  for  more  detail).  The  nonlinear  vector  function 
which  transforms  the  state  in  modified-polar  coordinates  to  the  state  in  Cartesian  coordinates  is  denoted 

§Mp(  — >R6  andisgivenby 
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Likewise,  the  nonlinear  vector  function  which  transforms  the  state  in  Cartesian  coordinates  to  the  state  in  modified- 
polar  coordinates  is  denoted  g£P (xc) :  R6  — >  R 8  and  is  given  by 
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The  above  equations  exactly  transform  the  state  estimate  from  either  coordinate  system  into  the  other  system.  What 
remains,  then,  is  the  transformation  of  the  covariance  matrices  from  one  coordinate  system  to  the  other.  It  is  this  step 


which  introduces  modeling  errors  and  is  only  approximate.  The  transformation  of  the  covariance  matrix  in  Cartesian 
coordinates,  Pc,  to  that  in  modified-polar  coordinates,  PMP  is  given  by  PMP  =  3qFPc^CF)T  where  J'PP  is  the 
Jacobian  matrix  of  the  transformation  from  Cartesian  to  modified-polar  coordinates  and  is  given  by 
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.  Finally,  the  Jacobian  of  the 


transformation  from  modified-polar  to  Cartesian  coordinates  is  given  by 
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3.3  Multiple  Hypothesis  Hybrid  Coordinate  Extended  Kalman  Filter  (MH-HC-EKF) 

A  common  problem  with  applying  the  EKF  to  the  bearings-only  tracking  problem  is  divergence  of  the  filter.  It  is 
observed  that  this  generally  occurs  with  poor  initialization.  The  range-parameterized  algorithm  simply  runs  several 
EKFs  in  parallel,  each  initialized  differently  [12,  13].  Typically,  the  jammer  is  known  to  lie  in  a  certain  region  of  the 
state  space  through  prior  information.  The  space  in  the  direction  of  the  LOS  is  chopped  up  into  range  bins  and  each 
filter  is  initialized  using  the  mean  of  each  range  bin  as  the  range  to  the  target.  Given  the  first  bearing  measurement,  the 
initial  estimate  of  the  jammer  position  for  each  EKF  is  determined  by  this  bearing  and  the  mean  value  of  the  range  bin 
for  the  corresponding  EKF.  Thus,  multiple  filters  are  estimating  the  position  of  the  jammer,  each  initialized  with  a 
different  value.  Each  filter’s  estimate  can  be  considered  as  a  hypothesis.  The  “probability”  that  each  filter’s  estimate  is 
the  correct  one  is  found  using  the  value  of  the  measurement  likelihood.  Note  that  this  is  somewhat  ad-hoc  in  the  sense 
that  the  hypotheses  are  not  disjoint  events.  Another  perspective  on  this  approach  is  that  we  have  several  independent 
filters  tracking  the  same  object.  The  measurement  likelihood  is  an  indication  of  how  well  the  filter’s  predicted 
measurement  matches  the  actual  measurement.  Thus,  the  measurement  likelihood  is  a  confidence  value  on  the  accuracy 
of  each  filter’s  estimate.  This  confidence  value  is  used  to  weight  each  filter  estimate.  The  overall  estimate  is  then  the 
weighted  sum  of  each  filter  estimate.  The  benefit  over  a  single  Kalman  filter  is  that  any  of  the  filters  which  begin  to 
diverge  will  exhibit  this  divergence  in  the  measurement  likelihood  leading  to  low  confidence  values.  Thus,  these  filters 
will  contribute  little  to  the  overall  estimate  and,  when  subject  to  a  pruning  technique,  will  be  discontinued. 

There  is  certainly  intuitive  appeal  to  this  multiple  EKF  approach.  We  implemented  this  for  the  hybrid  algorithm  but 
found  the  results  to  be  worse  than  for  a  single  hybrid  filter.  The  reasoning  is  that  the  hybrid  filter  rarely  diverges  (we 
found  no  cases  for  which  it  did)  and,  thus,  all  parallel  filters  continue  to  track  and  converge  to  the  correct  estimate. 
However,  depending  on  the  initial  estimate,  the  filters  have  different  convergence  rates,  and  the  good  estimates  get 
diluted  with  the  poor  estimates.  Perhaps  for  a  standard  EKF  which  is  prone  to  divergence  under  poor  initialization  the 
range-parameterized  algorithm  will  improve  performance. 


3.4  Particle  Filter  (PF) 

The  Bayesian  framework  for  filtering  problems  represents  a 
theoretically  sound  approach  which  entails  determining  the  PDF  of  an 
unknown  parameter.  The  following  formula,  referred  to  as  Bayes’ 
Rule,  provides  a  means  by  which  the  prior  PDF  of  some  unknown 
random  time -varying  parameter  x(k)  is  updated  using  a  likelihood 
function 


p(x(k)\Yk)  = 


p(y(k)\x(k))p(x(k)\Yk-i) 
p(y(k)  |  Yk~l) 

(BU) 


where  Yk  =  {y(l),y(2),...,y(£)}  represents  a  sequence  of 
measurements.  The  measurements  are  related  to  the  unknown 
parameter  through  the  following  measurement  equation 
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Figure  6:  HC-EKF  vs.  PLLS  with  3-DOF 
model  and  jammer  at  (-1,  0,  0)  km 


y(£)  =  MxW,v(£)) 


(NLME) 


where  h (■,■)  is  a  possibly  time -varying  nonlinear  function  of  the  unknown  parameter  x(/{) .  and  \(k)  is  typically 

assumed  to  be  a  sample  of  a  white  noise  measurement  process  whose  PDF  is  known.  The  state  is  assumed  to  evolve 
according  to  the  following  general  dynamic  equation 

x(k  + 1)  =  (x(£),  w(£))  (NLSE) 


where  //,  (•,•)  may  be  a  time-varying  nonlinear  function  of  the  state  and  w (k)  is  typically  modeled  as  a  white  noise 
process  with  a  known  PDF.  The  difficulty  in  analytically  determining  the  PDF  of  x(k)  is  that  only  for  special  cases  does 

a  closed-form  (i.e.,  finite)  description  of  p{x{k)  \  Yk)  exist.  When  (NLME)  and  (NLSE)  take  on  the  special  linear 
forms 


x(k  + 1)  =  F(£)x(£)  +  r(£)w(£) 
y(k)  =  H(£)x(&)  +  \(k) 


(LE) 


where  \\(k)  and  \(k)  are  samples  from  a  Gaussian  white  noise  process,  then  the  PDFs  are  all  Gaussian  and,  thus,  are 
completely  characterized  by  the  first  and  second  moments  (i.e.,  mean  and  covariance).  For  the  linear  case  modeled  as  in 
(LE),  the  resulting  Bayes’  implementation  is  the  well-known  Kalman  Filter.  Flowever,  for  a  plethora  of  other  cases 
(e.g.,  nonlinear,  non-Gaussian),  the  Kalman  filter  is  inappropriate.  While  some  success  has  been  achieved  by  applying 
the  EKF  (Extended  Kalman  filter),  derived  by  linearizing  the  measurement  equation  and/or  dynamic  equation,  the 
behavior  of  the  EKF  is  often  unpredictable.  The  difficulty  in  applying  the  EKF  lies  in  the  underlying  assumption  that 
the  PDF  of  the  unknown  parameter  (i.e.,  target  state)  can  be  sufficiently  approximated  by  a  Gaussian  density  function,  a 
dubious  assumption,  particularly  for  the  case  of  multimodal  PDFs. 

When  linearization  of  the  measurement  model  and/or  dynamic  model  is  not  feasible  or  leads  to  unsatisfactory  results, 
another  approach  to  the  implementation  of  Bayes’  methods  is  the  approximation  of  the  true  PDF  of  x(k)  using  a  model 
density.  Various  approximations  include  a  piece-wise  constant  approximation,  point-mass  approximation  or  Gaussian 
sum  approximation  of  the  true  posterior  PDF  of  the  target  state.  These  “deterministic”  methods  all  involve 
approximating  the  PDF  of  the  state  on  a  fixed  grid  over  the  state  space  and,  therefore,  suffer  from  the  “curse  of 
dimensionality.”  To  cope  with  the  computational  complexity  of  these  fixed-grid  approaches,  sub-optimal  pruning  of 
basis  functions  (i.e.,  points  in  the  grid)  is  usually  required.  In  contrast,  the  method  described  next  naturally  implements  a 
randomly  evolving  grid  based  on  the  likelihoods  of  the  received  data. 

Stochastic  sampling  methods  randomly  sample  points  in  the  state  space  and  calculate  the  PDF  at  these  points.  Referred 
to  generally  as  Monte  Carlo  methods,  the  computation  required  for  their  implementation  in  real-time  systems  precluded 


them  from  use  except  in  statistical  batch  processing  type  applications. 
With  the  revolution  in  computing  power  over  the  last  decade  and 
advent  of  parallel  computing  structures  (particle  filtering  methods  are 
easily  parallelizable),  Monte  Carlo  methods  have  experienced  a 
resurgence  in  the  field  of  estimation  and  filtering  [14].  Particle  filters 
(also  referred  to  as  weighted  bootstrap  filters)  are  based  on  the  solution 
to  the  following  problem,  “Given  samples  from  a  PDF  g(x ) ,  can  one 
generate  samples  of  an  arbitrary  PDF,  / (x)  ?”  There  are  different 
methods  by  which  this  can  be  achieved,  one  of  which  is  the  weighted 
bootstrap  method  on  which  particle  methods  are  based.  For  illustrative 
purposes  we  demonstrate  the  implementation  of  (BU),  which  is  the 
“update”  step  in  the  Bayes’  approach  to  filtering.  Given  that  we  have 
samples  xt(k) ,  i  =  l,2,...,n,  which  were  drawn  from  the  predicted 

density,  p(x(k)  \  Yk  '),  the  following  procedure  is  used  to  simulate  a 

set  of  samples  drawn  from  p{x{k )  Y1' ) .  First,  after  receiving  a 
measurement,  y  (k) ,  the  weights 
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Figure  7:  HC-EKF  vs.  PLLS  with  3-DOF 
model  and  jammer  at  (-2,  0,  0)  km 


p(y(k)  I  Xj(k)) 
Zp(y(k)\xj(k)) 


i  =  1,2,3,...,  n. 


are  calculated  using  the  likelihood  function  evaluated  at  y (k) .  Next,  a 
discrete  distribution  is  formed  over  the  samples  x;(£)with  qt  as  the 
weights.  Finally,  this  discrete  distribution  is  resampled  n  times  yielding 
a  set  of  points  which  simulate  samples  drawn  from  the  density, 

p(x(k)  Yk ) .  The  requirements  for  implementing  the  weighted 
bootstrap  filtering  algorithm  are  that  a  prior  distribution  on  the 
unknown  random  parameter  is  available  for  sampling,  the  likelihoods 
are  available  and  the  PDF  of  w (k)  in  (LE)  is  available  for  sampling. 
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Note  the  absence  of  reliance  on  linearity  or  Gaussian  noise  statistics.  Figure  8:  FIC-EKF  localization  error  with 
In  fact,  a  closed-form  expression  for  the  likelihood  function  is  not  even  6-DOF  model  and  moving  jammer 
necessary.  It  is  sufficient  that  the  likelihood  be  evaluated  pointwise  for 

any  point  in  the  state  space.  Thus,  particle  filtering  methods  are  applicable  to  a  much  wider  range  of  estimation 
problems  than  traditional  linear  Kalman  Filters.  Furthermore,  there  exists  a  justification  [14]  of  the  convergence  of 
particle  approximations  to  the  true  PDFs,  results  which  are  lacking  for  ad-hoc  methods  such  as  the  EKF. 


4.  SIMULATION  RESULTS 

All  four  algorithms,  the  PLLS  algorithm,  the  HC-EKF,  the  MH-HC- 
EKF,  and  the  PF,  were  individually  optimized  for  the  simulation 
assumptions  including  measurement  and  process  noise  values.  Each 
algorithm  was  run  for  varying  flight  times  and  conditions  and  with  both 
3-DOF  and  6-DOF  dynamic  models.  The  scenario  with  the  3-DOF 
munition  model  was  set  up  such  that  the  munition  lost  GPS  lock  at  (0,  - 
8,  7.62)  km  while  traveling  due  north  at  Mach  0.8.  The  unaided  flight 
time  is  approximately  40  seconds.  The  6-DOF  scenario  was  set  up  such 
that  the  munition  lost  GPS  lock  at  (0,  -15,  7.62)  km  while  flying  due 
north  at  Mach  0.8.  The  unaided  flight  time  was  approximately  73 
seconds.  We  found  all  four  algorithms  to  be  stable  but  we  did  verify  the 
bias  properties  the  PLLS  filter.  Essentially,  the  PLLS  localization  error 
was  highly  dependent  on  the  system  noise,  the  distance  to  the  jammer, 
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Figure  9:  HC-EKF  performance  for  3-DOF 
and  6-DOF  dynamic  models 


Figure  10:  MH-HC-EKF  for  various 
numbers  of  parallel  filters 


and  the  geometry.  The  performance  of  the  PLLS  algorithm  was  also 

degraded  when  the  jammer  was  moving  with  a  constant  velocity.  We  —  h=i 

now  present  results  that  are  representative  of  the  performance  of  each 
algorithm. 

Figure  6  shows  the  localization  performance  of  the  PLLS  algorithm 
and  the  HC-EKF  when  the  3 -DOF  munition  model  is  used  and  the 
stationary  jammer  is  at  (-1,  0,  0)  km.  We  note  that  localization  error  in 
meters  CEP  for  the  HC-EKF  was  less  than  the  error  for  the  PLLS 
algorithm  for  the  duration  of  the  flight  and  is  significantly  lower  at  the 
end  of  the  Eight  (<  100  m  CEP  vs.>  250  m  CEP).  The  difference  in  the 
performance  of  the  two  algorithms  was  significantly  more  pronounced 
when  the  jammer  was  moved  further  away  from  the  origin  (the  final 
destination  of  the  munition.)  Figure  7  demonstrates  this  fact  for  the 
same  initial  conditions  as  the  previous  figure  except  that  the  jammer  is 
now  at  (-2,  0,  0)  km.  The  final  localization  error  for  the  HC-EKF  is 
~  1 00  m  CEP  while  the  localization  error  for  the  PLLS  algorithm  is  greater  than  400  m  CEP. 

Figure  8  shows  the  response  of  the  HC-EKF  when  the  jammer  is  moving  with  a  constant  velocity  of  (v/V 2,  v/V2,  o) 
where  v  =  5,  10,  15,  and  20  m/s.  In  this  case  the  6-DOF  munition  is 
utilized.  The  figure  demonstrates  that  the  final  localization  error  for  the 
HC-EKF  is  less  than  100  m  CEP  for  all  cases.  Although  we  have  not 
included  supporting  data,  the  localization  error  for  the  HC-EKF  also 
did  not  appear  to  depend  on  the  range  to  the  jammer,  i.e.,  the  HC-EKF 
appears  to  be  unbiased  with  respect  to  range.  Figure  9  demonstrates 
that  the  filter  performance  did  not  significantly  depend  on  the  munition 
model  employed.  It  shows  the  HC-EKF  response  for  a  stationary 
jammer  located  at  (  -r,  0,  0)  for  r  =  1,  2,  and  4  km.  In  this  case  the  3- 
DOF  munition  had  the  same  initial  conditions  as  the  6-DOF  model 
with  a  73  second  unaided  flight  time. 

Figures  10  and  11  compare  the  performance  of  the  HC-EKF  to  that  of 
the  MH-EC-EKF  and  the  PF.  In  both  cases  the  3-DOF  munition  model 
is  employed.  In  Figure  10  we  varied  the  number  of  initial  hypothesis 
considered  for  the  range-parameterized  (multiple-hypothesis)  HC-EKF. 

The  figure  shows  the  response  of  the  algorithm  for  H  =  1,6,  11,  16,  21, 

and  26  where  H  is  the  number  of  parallel  filters  considered.  It  appears  that  a  single  HC-EKF  (H  =  1)  is  sufficient  for 
good  performance.  Figure  11  compares  the  performance  of  the  HC-EKF  with  that  of  the  PF.  The  figure  demonstrates 
that  while  the  final  localization  error  of  the  PF  is  comparable  under 
various  different  geometries  to  that  of  the  HC-EKF  (~  100  m  CEP),  the 
HC-EKF  is  able  to  maintain  a  significantly  lower  localization  error  for 
a  much  greater  portion  of  the  trajectory. 

Finally,  Figure  12,  shows  typical  impact  errors  for  the  case  when  the 
munition  pursued  the  jammer  instead  of  the  intended  target  at  the 
origin.  While  the  impact  errors  of  Figure  12  were  less  than  1.7  m  CEP, 
we  found  that  the  impact  errors  for  various  stationary  jammer 
geometries  and  various  unaided  flight  time  durations  were  all  within 
2.5  m  CEP  for  both  the  3-DOF  and  6-DOF  munition  models.  The 
guidance  laws  for  both  the  3-DOF  and  6-DOF  munition  models  were 
based  on  classical  proportional  guidance.  A  commanded  impact  angle 
of  71°  was  always  achieved.  Only  the  PLLS  and  HC-EKF  algorithms 
were  utilized  for  home-on-jam  scenarios.  Both  algorithms  performed 
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Figure  12:  Typical  impact  errors  after 
homing  on  the  jammer 


similarly  for  stationary  jammers  while  the  HC-EKF  performed  significantly  better  for  moving  jammers. 


5.  CONCLUSIONS 

We  tested  several  bearings-only  localization  algorithms  for  the  purpose  of  locating  GPS  jammers  from  a  munition  that 
has  lost  GPS  lock.  We  found  that  the  HC-EKF  performed  very  well  under  all  scenarios  that  we  considered.  While  the 
PLLS  filter  was  stable,  it  was  found  to  be  biased  as  a  function  of  range,  geometry  and  noise  level.  The  multiple- 
hypothesis  EKF  in  hybrid-coordinates  did  not  seem  to  provide  any  advantages  over  the  single  HC-EKF.  While  the  PF 
errors  all  did  converge  to  approximately  the  same  values  as  the  HC-EKF,  the  performance  of  the  HC-EKF  algorithms 
was  preferable  over  the  entire  munition  trajectory  and  because  it  is  computationally  less  demanding  than  the  PF. 
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